Electron correlation and the phase diagram of Si 
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Previous first-principles calculations of the melting properties of Si, based on the local-density 
approximation (LDA) for electronic exchange-correlation energy, under-predict the melting tem- 
perature by ~ 20 %. We present new first-principles results demonstrating that this problem is 
due to non-cancellation of exchange-correlation errors between the semiconducting solid and the 
metallic liquid. It is shown that other sources of error, particularly those due to system size and 
Brillouin-zone sampling, can be made negligible. The same LDA errors cause an underprediction 
of the pressure of the diamond-Si — > beta-tin-Si transition. The generalized-gradient approximation 
largely corrects both features of the Si phase diagram. 
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The long-standing ambition of calculating phase dia- 
grams from first-principles quantum mechanics has be- 
come a reality in the last 10 years G|— 0] . An impor- 
tant stimulus to the recent developments was the pa- 
per of Sugino and Car (hereafter SG) on the melting of 
Si The authors showed how the technique of ther- 
modynamic integration [0] combined with first-principles 
molecular dynamics (FPMD) based on density func- 
tional theory (DFT) can be used to calculated the free en- 
ergy of solids and liquids, and hence melting curves, with 
no experimental input apart from fundamental constants. 
But although their paper was influential, their numerical 
results on Si were not very satisfactory, since their pre- 
dicted melting temperature (T m — 1350 K) was ~ 20 % 
below the experimental value (1685 K) ||]. Our purpose 
here is to identify the cause of this discrepancy, which we 
shall argue comes from non-cancellation of DFT errors 
between the solid and liquid phases, and specifically from 
errors of the local-density approximation (LDA) used by 
SG. This has implications for the reliability of other first- 
principles work on phase diagrams. 

The basic approximation in any DFT calculation is the 
algorithm adopted for exchange- correlation energy E xc . 
Provided one can eliminate all other sources of error in 
calculating total energies and doing the statistical me- 
chanics, then failure to reproduce experimental melting 
properties must be due to errors in E xc . But it is often 
claimed, even by first-principles practitioners, that these 
other sources of error cannot be made small enough; in 
particular, it is claimed that first-principles calculations 
cannot yet be performed on large enough systems to ren- 
der size errors negligible Q. This was one of the major 
issues addressed by SG, who made strenuous efforts to en- 
sure that their non-£" xc errors were negligible; the results 
presented later indicate that they were largely success- 
ful. Turning to E xc errors, the crucial question is the ex- 
tent to which they cancel between the coexisting phases. 
Since diamond-structure Si (d-Si) is a four-fold coordi- 
nated semiconductor and liquid Si (1-Si) is an approxi- 
mately six- fold coordinated metal B], electron screening 



is likely to be very different in the two phases, so that 
non-cancellation of E xc errors becomes an important is- 
sue. In considering this, we arc helped by the fact that 
the pressure-stabilized /3-tin structure (/3-tin-Si) closely 
resembles the liquid in being metallic and six-fold coordi- 
nated. This suggests that there should be a close relation 
between the effect of E xc errors on the melting tempera- 
ture and on the d-Si — > /3-tin-Si transition pressure, and 
an analysis of this relation will help us to confirm that 
errors in the LDA representation of E xc account for the 
under-prediction of T m . 

Our first-principles calculations employ Vanderbilt 



ultra-soft pseudopotentials |10| and plane- wave basis 
sets. Most of our calculations are based on the local- 
density approximation (LDA) for E xc used by SG, but we 
shall also present results using the generalized-gradient 
approximation (GGA) The calculations were done 

with the VASP code [ Ejf . The plane-wave cut-off was 
150 eV, which gives a convergence of 6 meV/atom in 
the difference of total (free) energies between liquid and 
solid, and the pseudopotential core radii were 1.31 A. 
Our strategy for computing the free energies of solid and 
liquid differs somewhat from that of SG, and closely re- 
sembles that used in our recent work on Fe and Al || . 

The Hclmholtz free energy F of the solid can be writ- 
ten as F = Fp Cr f + F^ib , where F pel { is the free energy of 
the perfect non-vibrating crystal (it is a free energy, be- 
cause we allow for thermal electronic excitations), and 
FVib is the contribution from lattice vibrations. The 
latter is written as F v ib = Fh arm + F an h a rm- The har- 
monic free energy per atom Fharm in the classical limit 
(melting occurs well above the Debye temperature) is: 
-Fharm = Sk^T \t\(Tiuj / k-&T) , where the geometric-mean 
frequency Q is given by: 



ln(^HA^]TlnK s ), 



(1) 



k.s 



with the sum going over wavevectors k and branches s in 
the Brillouin zone, N^ s being the number of terms in the 
sum. The phonon frequencies lu^ are calculated using 
the small-displacement method. [Ol 
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The anharmonic contribution F an harm turns out to be 
very small (typically ~ 15 meV/atom near the melting 
temperature) , so that it is accurately given by the second- 
order expansion: 

F a nharm — (^anharm)harm (^anharm)harm/2fcBF ; (2) 

where ^7 an harm is the anharmonic part of the first- 
principles total energy, and the thermal averages ( • )harm 
are evaluated in the canonical ensemble of the first- 
principles harmonic system. We have verified the accu- 
racy of Eq. (^|) by comparing it with the exact expression: 

fanharm = ~k B T ln(exp(- [/ a iiharm/fcBF) ) harm . 

Our calculations of F pcr f were performed on the 
primitive two-atom unit cell, at volumes from 16 to 
22 A 3 /atom with fc-point sampling dense enough to give 
a precision of ~ 0.1 meV/atom. Results were fitted to the 
Birch-Murnaghan jl4} form, which reproduces the data 
to within ~ 0.1 meV/atom. For -Fhanni we calculated 
the force-constant matrix using 54-atom cells, with spot- 
checks on cells of up to 250 atoms indicating convergence 
to within ~ 2 meV/atom. The calculations were done at 
volumes from 18 to 21 A 3 /atom, and In (a;) was fitted to 
a second-order polynomial ln((D) = a + bV + cV 2 , which 
gives a fitting error in F^ aim of ~ 1 meV/atom. The ther- 
mal averages needed to calculate F an harm were done on 
a 54-atom cell at volumes of V = 18 and 20 A 3 /atom 
and temperatures of 1000, 1500 and 2000 K. The re- 
sults are accurately reproduced by the quadratic form 
-Fanharm = aT 2 , and to the accuracy we require it is 
enough to take the value a = 7 x 10~ 9 eV K~ 2 for both 
volumes. 

The free energy of the liquid is calculated using 
thermodynamic integration (TI), with the Stillinger- 
Weber |H| empirical total-energy model used as refer- 
ence system. The difference of Helmholtz free energy 
AF = Fai — F Te f between the ab initio and Stillinger- 
Weber systems is obtained using the standard formula: 

AF = f dX (U A i - U rci ) x , (3) 
Jo 

with Uai and U TO { the ab initio and reference total-energy 
functions, and ( • }\ the thermal average evaluated in the 
ensemble of the system whose total-energy function is 
U\ = (1 — X)U Te f + XUai- In practice, the integral over A 
is performed either by evaluating (U — U TC {)\ at a set of 
A values and using Simpson's rule, or by using 'adiabatic 
switching', in which A is slowly and continuously varied 
between the two limits ]16|] . The reference free energy 
F Te f is calculated by thermodynamic integration starting 
from the Lennard- Jones system, for which accurate free 
energies have been published |l7| . The calculation of F re f 
is done on very large systems, so that it is converged with 
respect to system size to better than 1 meV/atom. 

As has often been stressed M , the final results for Fai 
do not depend on the choice of reference system, but the 



efficiency of the calculations can be greately improved 
by careful tuning of the reference system, the criterion 
being that the strength of the fluctuations of Uai — U le { 
should be made as small as possible. We find that in 
this sense the original parameters of the Stillinger- Weber 
model are far from optimal for liquid Si. Using ab initio 
MD simulations of 1-Si at the state V — 18.16 A 3 /atom, 
T = 2000 K, we have varied the Stillinger- Weber pa- 
rameters to minimise the fluctuation strength, and it is 
the resulting 'optimized' SW model that we use as our 
reference model. 

We made thorough tests of the convergence of AF 
with respect to system size and electronic fc-point sam- 
pling by calculating it at the representative state point 
V = 17 A 3 /atom and T = 1750 K, using systems of up 
to 512 atoms and up to 36 Monkhorst-Pack |l8| fc-points 
(results of these tests in Table 1). The tests were done 
as follows. The T-point results were obtained by explicit 
simulations on systems of all sizes, with AF calculated 
by thermodynamic integration (Eq. @). In most cases, 
we used the five A values 0.0, 0.25, 0.5, 0.75 and 1.0 to- 
gether with Simpson's rule, and comparisons with other 
sets of A values show that the residual error from the in- 
tegration itself is less than 5 meV/atom. We then used 
thermodynamic integration, with the T-point system as 
the reference system, to obtain the results for other fc- 
point samplings. For systems of N > 128 atoms, the 
fluctuations of the difference of energies calculated with 
T-point and more fc-points is small enough to allow the 
second-order expansion to be used instead of explicit TI, 
but for N = 64 this is not adequate and we used explicit 
TI. 

The results of Table 1 show that with 64 atoms and 
four fc-points the free-energy difference AF between ab 
initio and optimized Stillinger- Weber is converged to bet- 
ter than 10 meV/atom, and we have used this system to 
obtain Fai for the liquid at the set of state points V = 16, 
17, 18, 19 and 20 A 3 /atom and T = 1250, 1500 and 
1750 K. At each T, Fai was fitted to a Birch-Murnaghan 
equation of state, the residual fitting error being no more 
than 2 meV/atom. 

Our fitted ab initio Helmholtz free energies of d-Si 
and 1-Si allow us to obtain the Gibbs free energy G = 
F — V(dF/dV)T, and hence the melting curve. The zero- 
pressure results for T m and the entropy and volume of 
fusion, AS* and AV, are compared in Table 2 with those 
of SG and the experimental values. Our very close agree- 
ment with the SG value of T m (difference of only 50 K) 
confirms that their size and fc-point errors were indeed 
very small, and also confirms that LDA under-predicts 
T m by ~ 20 %. We note that our AS* and AV values are 
both somewhat greater than those of SG. 

We now turn to the matter of non-cancelling LDA er- 
rors between phases, exploiting the electronic and struc- 
tural similarity between d-Si and /3-tin-Si. At room tem- 
perature, the transition d-Si — > /3-tin-Si occurs at an ex- 
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perimental pressure in the range 10.3 - 12.5 GPa |J] 
(although also a low value of 8.8 GPa has been re- 
ported ^0|). Earlier LDA calculations on the static 
zero-temperature crystals gave transition pressures in the 
range 7.8—8.4 GPa pl"| , p2"| , and our own calculations yield 
the value 7.8 GPa, which agrees closely with the earlier 
values. However, it has been shown that temperature 
has a strong influence on the transition pressure, which 
drops by ~ 20 % as T goes from K to room temper- 
ature |^3| , so that the temperature-corrected LDA pres- 
sure is too low by at least 4 GPa. It is also known that 
the generalized-gradient approximation (GGA) for E xc 
significantly improves the predicted transition pressure. 
With the Perdew-Wang GGA @, we find a transition 
pressure of 11.7 GPa (9.4 GPa when corrected to room 
temperature) , which agrees closely with earlier GGA val- 
ues f22| . Our calculations show that the main reason why 
LDA under-predicts the transition pressure is that it er- 
roneously shifts the energy of d-Si upwards relative to 
/3-tin-Si. The GGA goes a long way towards correcting 
this destabilization of d-Si. But a low melting tempera- 
ture is also a sign of an erroneous destabilization of d-Si, 
and we hypothesize that the same underlying E xc error 
is responsible for both under-predictions. 

To test this hypothesis, we have recalculated the melt- 
ing properties using GGA. It is instructive to do this by 
evaluating the free energy difference between the LDA 
and GGA systems. We have therefore performed long 
simulations for solid and liquid at the zero pressure vol- 
umes using the LDA, and calculated the GGA energies at 
a number of statistically independent configurations, for 
both the solid and the liquid. The calculations have been 
done on cells containing 64 atoms with four /c-points, 
and spot-checked with calculations on cells containing 
512 atoms and T-point sampling. Firstly, we found that 
the energy differences between GGA and LDA are basi- 
cally constant, i.e. do not depend on the configurations of 
the atoms, which confirms the idea that the shift should 
be the same as for the low temperature static lattices. 
Secondly, we found that the free energy of the liquid is 
raised by 88 meV/atom relative to that of d-Si. Given an 
LDA entropy change on melting of 3.5/cB/atom, it is easy 
to work out a shift of melting temperature GGA-LDA of 
292 K, bringing the GGA result to 1590 K, in much closer 
agreement with the experimental datum. We also found 
that, at the volumes corresponding to the LDA zero pres- 
sure, the GGA pressures are about 3.5 GPa larger than 
the LDA ones, so the GGA zero pressure volumes are 
larger. However, the bulk moduli for the solid and the 
liquid at the melting temperature are 78 and 34 GPa re- 
spectively, so the liquid will expand more than the solid 
in the GGA. We can estimate a new volume change on 
melting of 9.4%, which is also in somewhat better agree- 
ment with the experiments. 

In summary, we have shown that the key issue in a first- 
principles account of the melting properties of Si is non- 



cancellation of exchange-correlation errors between solid 
and liquid because of their different electronic structure. 
Technical errors due to system size and k- point sampling 
are readily brought under tight control. The basic reason 
why this can be done is that system size affects only the 
small difference of free energy between the first-principles 
system and a carefully designed reference system. The 
non-cancellation of exchange-correlation errors between 
coexisting semiconductor and metal is also responsible for 
difficulties in predicting the pressure of the diamond-Si — > 
/3-tin-Si transition, and there is a quantitative relation 
between the error in this transition pressure and the error 
in melting temperature. The general implication is that 
for phase equilibria in which the coexisting phases have 
essentially the same electronic structure, e.g. the melting 
of high- pressure Fe, DFT calculations can be expected 
to predict phase equilibria with satisfactory accuracy. 
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AF 32 


AF 36 


64 


-4.165(5) 


-4.262(5) 


-4.253(5) 


-4.257(5) 


-4.257(5) 


128 


-4.282(5) 


-4.250(5) 








216 


-4.281(5) 


-4.262(5) 








512 


-4.248(5) 


-4.251(5) 









TABLE I. Difference AF of Helmholtz free energy (in 
units of eV/atom) at thermodynamic state V = 17 A 3 /atom, 
T = 1750 K, between the ab-initio and the Stillinger- Weber 
potential as function of size of simulated system (number of 
atoms N) and number of Monkhorst-Pack k-points (subscript 
on AF). 



This work (LDA) This work (GGA) Sugino and Car || Expei 



T m (K) 

AVm/V, 

AS m 
dT m /dP 



1300(50) 
0.142 
3.5 
-58 



1590(50) 
0.094 



1350(100) 
0.1 
3.0 
-50 



a Rcf. pr 

b Ref. [24 1 
c Ref. [25 
d Ref. M 

TABLE II. Comparison of calculated and experimental 
melting properties of Si at ambient pressure: melting tem- 
perature T m , volume change AV m divided by volume of solid 
at melting temperature, entropy change AS m per atom di- 
vided by Boltzmann's constant, and slope of melting curve 
dT m /dP (units of K GPa" 1 ). 
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